# bar plot with custom error bars
# 23 February 2013
# NJG
opar <- par(no.readonly=TRUE)
Ant.Data <- read.csv("antcountydata.csv", header=TRUE)
attach(Ant.Data)

ecoregion.means <- by(n.species,ecoregion,mean)

barplot(ecoregion.means,names.arg=letters[1:5],
        ylim=c(0,100))
midpoints <- barplot(ecoregion.means,plot=FALSE)
ecoregion.sd <- by(n.species,ecoregion,sd)

arrows(midpoints,ecoregion.means,midpoints,
       ecoregion.means+ecoregion.sd,angle=90,
       length=0.1)

# For a double-headed arrow, adjust starting
# points and use code=3
arrows(midpoints,ecoregion.means-ecoregion.sd,midpoints,
       ecoregion.means+ecoregion.sd,angle=90,
       length=0.1,code=3)
points(midpoints,ecoregion.means,pch=19)

# plotting points with asymmetric
# confidence intervals
par(opar)

# using the quantile function
z <- rnorm(1000)
quantile(z,probs=0.95)
hist(z)
abline(v=quantile(z,probs=c(0.05,0.95)),lty="dashed",col="red")
abline(v=quantile(z,probs=c(0.025,0.975)),lty="dotted",col="blue")

# Calculate Confidence intervals with the by function
quantile(n.species,probs=0.025)
Con.low <- by(n.species,ecoregion,quantile,probs=0.025)
Con.high <- by(n.species,ecoregion,quantile,probs=0.975)

# set up a blank plot (show code first)
plot(1:5,1:5)
plot(1:5,1:5,xlim=c(0.5,5),ylim=c(0,100))
plot(1:5,1:5,xlim=c(0.5,5),ylim=c(0,100),type="n")

# fully blank plot
plot(1:5,1:5,xlim=c(0.5,5),ylim=c(0,100),type="n",
     ,axes= FALSE,ann=FALSE)

# add axes with labels
axis(1,at=seq(1:5),labels=letters[1:5])
axis(2)
mtext("Number of Species",side=2,line=2.5)
box(bty="L")


arrows(1:5,ecoregion.means,1:5,Con.low,angle=90,length=0.1)
arrows(1:5,ecoregion.means,1:5,Con.high,angle=90,length=0.1)
points(1:5,ecoregion.means,pch=23,cex=1.5,bg="gold")
abline(h=seq(0,100,by=20),col="gray",lty="dotted")

# add minor tick marks
# first find your location in the plot
points(1,0)
points(0,0)
points(0.5,0)
par("usr") # get limits

# try plotting a point points(0.32,10)
points(0.32,10)
par("xpd")
par(xpd=TRUE) # change clippig to figure region
points(0.32,10)
# first refresh plot
segments(rep(0.32,5),seq(10,90,length=5),rep(0,5),seq(10,90,length=5))
  
############################################
detach(Ant.Data)
par(opar)
rm(list = ls())